/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
* PLikeliHood.java
*
* Created on November 22, 2007, 9:11 AM
*
* To change this template, choose Tools | Template Manager
* and open the template in the editor.
 */
package genomemap.worker;

import genomemap.data.DataSource;
import genomemap.data.PDataSource;
import genomemap.worker.helper.PhysicalMapHelper;

import gui.genomemap.inputModel.ComPLHoodIModel;

import javautil.collections.ArrayUtil;

/**
 * A class dedicated to the physical likelihood of
 *  probe/clone orders. Note: This is a negative-log likelihood.
 */
public class PLikeliHood extends LHood {

    /**
     * physical data source
     */
    private PDataSource pData = null;

    /** Field description */
    private double[] probeSpacing = null;

    /**
     * Constructor
     */
    public PLikeliHood() {}

    /**
     * Constructs ...
     *
     * @param inputModel description
     */
    public PLikeliHood(ComPLHoodIModel inputModel) {
        setPData(DataSource.getPDataSource(inputModel.getLinkageGroup(), inputModel.getOrder()));
    }

    public static void main(String[] args) {

        PLikeliHood pHood = new PLikeliHood();

        pHood.setPData(DataSource.getPDataSource(1, new int[] { 1, 2, 3, 4, 5 }));
        pHood.setVerbose(true);
        pHood.runPHood();
    }

    public double runPHood() {

        setStatusInfo("Physical Likleihood Computation has Started");
        setProgress(10);

        int probe_count = pData.getOrder().length;
        int clone_count = pData.getCloneCount();

        /**
         * data is re-ordered according to the order in pData.
         */
        int[][] data           = pData.getdata();
        double probe_false_pos = pData.getProbFalsePos();
        double probe_false_neg = pData.getProbFalseNeg();
        int ch_length          = pData.getChLength();
        int clone_length       = pData.getCloneLength();

        /**
         * Note the data has been re-ordered according to the order in pData
         * the following probeOrder is just a dummy sequence from 1 to probeCount
         */
        int[] probeOrder = ArrayUtil.IntegerSequence(probe_count, false);

        /**
         * joint[i][j] = # of clones that hybridize to both i and j,
         * and joint[i][i] = 0, for all i
         * DEPENDENT PARAMS: probe_count,clone_count,data
         */
        int[][] joint = new int[probe_count][probe_count];

        for (int probeIndex1 = 0; probeIndex1 < probe_count; probeIndex1++) {

            for (int probeIndex2 = probeIndex1 + 1; probeIndex2 < probe_count; probeIndex2++) {

                for (int cloneIndex = 0; cloneIndex < clone_count; cloneIndex++) {

                    if ((data[cloneIndex][probeIndex1] == 1)
                            && (data[cloneIndex][probeIndex2] == 1)) {
                        joint[probeIndex1][probeIndex2]++;
                    }
                }

                joint[probeIndex2][probeIndex1] = joint[probeIndex1][probeIndex2];
            }
        }

        /**
         * aa is represented by matrix a(i,j) in the paper.
         * DEPENDENT PARAMS: probe_count,clone_count,data,probe_false_pos,probe_false_neg
         */
        double[][] aa = new double[clone_count][probe_count + 1];

        for (int cloneIndex = 0; cloneIndex < clone_count; cloneIndex++) {

            for (int probeIndex = 0; probeIndex < probe_count + 1; probeIndex++) {

                if (probeIndex == 0) {
                    aa[cloneIndex][probeIndex] = 0.0;
                } else if (data[cloneIndex][probeIndex - 1] == 0) {
                    aa[cloneIndex][probeIndex] = probe_false_neg / (1 - probe_false_pos);
                } else {
                    aa[cloneIndex][probeIndex] = (1 - probe_false_neg) / probe_false_pos;
                }
            }
        }

        /**
         * P = total # of 1`s in data: used in calculation
         * DEPENDENT PARAMS: probe_count, clone_count, data
         */
        int P = 0;

        for (int cloneIndex = 0; cloneIndex < clone_count; cloneIndex++) {

            for (int probeIndex = 0; probeIndex < probe_count; probeIndex++) {

                if (data[cloneIndex][probeIndex] == 1) {
                    P++;
                }
            }
        }

        /**
         * gap = total gap not filled by the probes; const1: used in the calculation
         * DEPENDENT PARAMS: ch_length, clone_length, probe_count, clone_count, probe_false_pos
         */
        double const1 = clone_count * Math.log(ch_length - clone_length)
                        - P * Math.log(probe_false_pos / (1 - probe_false_pos))
                        - probe_count * clone_count * Math.log(1 - probe_false_pos);
        int gap = ch_length - (probe_count * clone_length);

        /**
         * y = vector denoting spacings between the probes including
         * the left-most and right-most.
         * DEPENDENT PARAMS:
         */
        setStatusInfo("Computing initial probe spacings");
        setProgress(30);

        double[] y = PhysicalMapHelper.Initial_Ys(probeOrder, probe_count, clone_count, ch_length,
                         clone_length, gap, joint);

        /**
         * R = vector of constants used in the calculation
         * DEPENDENT PARAMS: ch_length,clone_length,probe_count,clone_count,aa,probeOrder
         */
        double[] R = PhysicalMapHelper.getR(ch_length, clone_length, probe_count, clone_count, aa,
                         probeOrder);

        /**
         * Conj(): STEEPEST DESCENT ALGORITHM
         *  GIVEN A VECTOR Y, MOVE TO THE MINIMUM ALONG THE GRADIENT DIRECTION
         * DEPENDENT PARAMS:
         */
        setStatusInfo("Running Steepest Descent Algorithm");
        setProgress(70);
        setLhood(PhysicalMapHelper.Conj(y, probe_count, clone_count, aa, probeOrder, clone_length,
                                        const1, R));
        setProbeSpacing(y);

        if (isVerbose()) {
            System.out.println("The log-likelihood is " + getLhood());
        }

        setStatusInfo("The Physical Log-likelihood is " + getLhood());
        setProgress(100);

        return getLhood();
    }

    /**
     * computes the function which is used for simulated annealing
     * Since Plikelihood by default returns the negative-loglikelihood
     * it is not futher negativitized. This will used as part of Annealable
     * interface to Simulated Annealing to serach for an order with minimum value
     * and thus max. loh-likelihood.
     */
    @Override
    public double compute() {

        getPData().setOrder(getOrder());

        return runPHood();
    }

    @Override
    public double initialValue() {
        return compute();
    }

    /**
     * physical data source
     */
    PDataSource getPData() {
        return pData;
    }

    /**
     * physical data source
     *
     * @param newVal
     */
    @Override
    public void setPData(PDataSource pData) {

        this.pData = pData;

        setOrder(pData.getOrder());
    }

    public double[] getProbeSpacing() {
        return probeSpacing;
    }

    void setProbeSpacing(double[] probeSpacing) {
        this.probeSpacing = probeSpacing;
    }

    @Override
    public double getLogLikeliHood() {
        return -getLhood();
    }
}    // end of class

